# cleanup names
library(janitor)
# tidyverse for general data manipulation
library(tidyverse)
# use to create time series tsibbles for fable package
library(tsibble)
# feasts includes forecast accuracy functions
library(feasts)
# easier manipulation of date/time data
library(lubridate)
# pretty display of plots
library(patchwork)
# mapping functions
library(purrr)
# rjhyndmans forecasting package
library(fable)
# better view of result tables
library(knitr)
# # last three are to create network diagrams for product hierarchy
library(treemap)
# more network diagram
library(data.tree)
# same as above
library(networkD3)
# more forecasting tools
library(fabletools)
# pretty output tables
library(kableExtra)Sample Hierarchical Forecasting
1 Summary
Using a sample data set from Kaggle, I’ll be analyzing sales history and developing forecasts at appropriate levels of aggregation. Models and forecasts will be generated using the Fable package from Rob J Hyndman detailed in Forecasting: Principals and Practice 3rd ed
I’ll be demonstrating the use of hierarchical forecasts against the sample data. Many businesses use a bottoms up approach, forecasting at the lowest level of aggregation (typically SKU), and summing up those forecasts to get totals for regions, categories, etc. I’ve had experience and frustration with this method, going through bottom level forecasts one by one, only to total them up to a number that doesn’t quite make sense given total/category level history (too low, too high, missing seasonality).
The sample data I selected was built for practicing analysis in Python, specifically for profit margin, sales metrics, and etc. I’ve found that sample data for forecasting tends to be very forecastable (low variability, few products, clear seasonality/trend), while actual business data tends to be much messier. I wanted to use a data set that was a little more painful to analyse that better represented challenges in business forecasting.
This workbook ended up being quite long since I built it while learning and working through the forecasting problem, rather than creating this doc as a showcase/presentation. I did my best to keep things organized and explain the process as I went though.
2 Setup
First step is to load the appropriate R packages, load our data, and get our first view of what we’re working with.
2.1 Load packages
Load R packages, code below includes quick description for each.
2.2 Load Data
I simply downloaded it as a .xls file and uploaded it from my PC to keep things simple for myself, so the process is just reading the excel file from my documents.
I loaded the file and then used the janitor package to change the column names to snake_case (my preference). First view of our data shows we have 9994 records across 21 columns. For forecasting purposes only a few of these will be useful, but we’ll keep them in for now until we can take a look at them.
There are two date type columns available, order_date and ship_date. I chose to use order_date for forecasting because it is closest to customer demand; ship_date can have many extra variables based on stock availability, warehouse issues, etc.
data <- readxl::read_xls("C:/Users/andre/Documents/r4ds/SampleSuperStore/Sample - Superstore.xls")
# clean up column names, using snake_case as default
data <- data |> janitor::clean_names()
# check column names
glimpse(data)Rows: 9,994
Columns: 21
$ row_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
$ order_id <chr> "CA-2016-152156", "CA-2016-152156", "CA-2016-138688", "U…
$ order_date <dttm> 2016-11-08, 2016-11-08, 2016-06-12, 2015-10-11, 2015-10…
$ ship_date <dttm> 2016-11-11, 2016-11-11, 2016-06-16, 2015-10-18, 2015-10…
$ ship_mode <chr> "Second Class", "Second Class", "Second Class", "Standar…
$ customer_id <chr> "CG-12520", "CG-12520", "DV-13045", "SO-20335", "SO-2033…
$ customer_name <chr> "Claire Gute", "Claire Gute", "Darrin Van Huff", "Sean O…
$ segment <chr> "Consumer", "Consumer", "Corporate", "Consumer", "Consum…
$ country <chr> "United States", "United States", "United States", "Unit…
$ city <chr> "Henderson", "Henderson", "Los Angeles", "Fort Lauderdal…
$ state <chr> "Kentucky", "Kentucky", "California", "Florida", "Florid…
$ postal_code <dbl> 42420, 42420, 90036, 33311, 33311, 90032, 90032, 90032, …
$ region <chr> "South", "South", "West", "South", "South", "West", "Wes…
$ product_id <chr> "FUR-BO-10001798", "FUR-CH-10000454", "OFF-LA-10000240",…
$ category <chr> "Furniture", "Furniture", "Office Supplies", "Furniture"…
$ sub_category <chr> "Bookcases", "Chairs", "Labels", "Tables", "Storage", "F…
$ product_name <chr> "Bush Somerset Collection Bookcase", "Hon Deluxe Fabric …
$ sales <dbl> 261.9600, 731.9400, 14.6200, 957.5775, 22.3680, 48.8600,…
$ quantity <dbl> 2, 3, 2, 5, 2, 7, 4, 6, 3, 5, 9, 4, 3, 3, 5, 3, 6, 2, 2,…
$ discount <dbl> 0.00, 0.00, 0.00, 0.45, 0.20, 0.00, 0.00, 0.20, 0.20, 0.…
$ profit <dbl> 41.9136, 219.5820, 6.8714, -383.0310, 2.5164, 14.1694, 1…Dates were stored as date times instead of just dates. We don’t need the time component, so we’ll change them to dates so they’ll easier to work with. This could have been addressed as the data was loaded from excel, but I find it easier to let the function do its work and clean up afterwards. I also dropped row_id, it’s not meaningful and just an index provided by the base data.
data <- data |>
mutate(order_date = as.Date(order_date), ship_date = as.Date(ship_date), .keep = "all") |>
select(-row_id)2.3 Structure of data
Using n_distinct as a function across all columns, we get a good first glance at the values, and can see which might be useful for creating meaningful slices. There are few enough ship modes for it to be useful, but the method of shipping is unlikely to have any forecasting value. Segment, region, category, and subcategory all have few enough identities for them to be useful so we’ll explore them.
I drop the following columns:
order_id, ship_date, country, ship_mode, customer_id, customer_name, country, city, state, postal_code, region, sales, discount, profit
Something is concerning with this view though, 1862 distinct values for product_id! That’s a lot of different items to forecast. What’s even more difficult is that there are only 14 distinct quantities listed, meaning that in our four years of daily sales history, there are only 14 different quantities sold per item. We’ll look at the quantities sold later, but this may indicate low volume/intermittent demand at the product_id level.
But our data is at the daily level, and is supposed to be “store sales” which should have many small purchases instead of larger quantities sold at once (like if we were forecasting shipments to distributors or vendors) and in the next step we’ll transform it to monthly and look at the distribution of quantities. (5009 order_ids is concerning compared to the 1334 ship dates, that’s only 3.8 shipments a day).
kable(
sapply(data, function(x) n_distinct(x)),
col.names = "distinct values"
)| distinct values | |
|---|---|
| order_id | 5009 |
| order_date | 1237 |
| ship_date | 1334 |
| ship_mode | 4 |
| customer_id | 793 |
| customer_name | 793 |
| segment | 3 |
| country | 1 |
| city | 531 |
| state | 49 |
| postal_code | 631 |
| region | 4 |
| product_id | 1862 |
| category | 3 |
| sub_category | 17 |
| product_name | 1850 |
| sales | 6144 |
| quantity | 14 |
| discount | 12 |
| profit | 7545 |
Also a quick check to see what our date range is, four years of demand, not too bad! This will let us use three years for testing, keeping 2017 as a holdout to test our forecasting methods
kable(range(data$order_date),
col.names = "",
caption = "min/max dates"
)| 2014-01-03 |
| 2017-12-30 |
3 Exploratory Data Analysis
Let’s get this data summarized by month, segment, category, sub_category, and product_id. I chose month because from my experience it’s the most common time period for business forecasts, and dropped region as an executive decision. If there was a business need, sure we could try to include it in the forecast models, but we’ll quickly see some evidence that slicing the data into smaller groups is likely unhelpful.
After the data is summarized into monthly buckets (instead of daily) we still only see 19 different quantities being sold, reinforcing the suspicion of intermittent demand at the product_id level.
We’ll come back to this issue, but lets check our other variables before we go deeper.
data <- data |>
group_by(
month = yearmonth(order_date),
segment,
product_id,
category,
sub_category
) |>
summarize(quantity = sum(quantity), .groups = "drop")
kable(
sapply(data, function(x) n_distinct(x)),
col.names = "distinct values"
)| distinct values | |
|---|---|
| month | 48 |
| segment | 3 |
| product_id | 1862 |
| category | 3 |
| sub_category | 17 |
| quantity | 19 |
3.1 Hierarchical structure
We need to check to make sure category/sub_category/product_id have strict nested relationship, meaning that each sub_category has a unique slice of the products, and each category has a unique slice of sub_categories.
By summarizing at product_id and sub_category we do get unique combinations, and filtering to product_ids and sub_categories with a count > 1 returns the empty tables below, proving a strict parent/child relationship.
I.E., product_id x is only in sub_category y, and sub_category y is only in category z. This is good!
A couple empty tables isn’t the best visualization for this, but it’s a good rigorous check.
data |>
group_by(category, sub_category, product_id) |>
summarize(n = n(), .groups = "drop") |>
group_by(product_id) |>
summarise(n = n(), .groups = "drop") |>
filter(n > 1)# A tibble: 0 × 2
# … with 2 variables: product_id <chr>, n <int>data |>
group_by(category, sub_category) |>
summarize(n = n(), .groups = "drop") |>
group_by(sub_category) |>
summarise(n = n(), .groups = "drop") |>
filter(n > 1)# A tibble: 0 × 2
# … with 2 variables: sub_category <chr>, n <int>Another way to check this and view the mapping is with a visual network. Feel free to mouse over, click and drag! For this visual I just wanted to practice making mappings and experimenting with the chart, but it’s not really more helpful than just a table.
product_heirarchy <- data |>
group_by(category, sub_category) |>
summarise(quantity = sum(quantity), .groups = "drop")
product_heirarchy$pathString <- paste("Total", product_heirarchy$category, product_heirarchy$sub_category, sep = "|")
asdf <- as.Node(product_heirarchy, pathDelimiter = "|")
asdf <- ToDataFrameNetwork(asdf, "name")
simpleNetwork(asdf, fontSize = 20)Unlike our strict product hierarchy however, segment is a grouping based on another variable, likely customer and their business purpose. Each segment represents sales from all 3 categories, and based off the unique product_id count for each segment, there is clearly some overlap (remember we have 1862 product_ids)
kable(
data |> group_by(segment, category) |>
summarize(n = n(), .groups = "drop") |>
group_by(segment) |>
summarise(category_count = n(), .groups = "drop"),
caption = "Number of categories"
)| segment | category_count |
|---|---|
| Consumer | 3 |
| Corporate | 3 |
| Home Office | 3 |
kable(
data |> group_by(segment, product_id) |>
summarize(n = n(), .groups = "drop") |>
group_by(segment) |>
summarise(product_id_count = n(), .groups = "drop"),
caption = "Number of products"
)| segment | product_id_count |
|---|---|
| Consumer | 1716 |
| Corporate | 1456 |
| Home Office | 1112 |
3.2 Shape of Demand
Now that we’ve got an understanding of the structure of the data, lets take a look at some more interesting charts and see how our demand is shaped.
3.2.1 Product ID
Using a histogram we can see the distribution of units sold per month per product, the most common values being 3. While this on its own would be relatively challenging to forecast, there’s another issue that stands out: there are no periods with zero sales! Our data is a record of sales, and when there is no sale, of course there is no record. To forecast however, we’ll need to fill in these zeros
monthly_prod_data <- data |>
group_by(month, product_id) |>
summarise(quantity = sum(quantity), .groups = "drop")
m <- round(median(monthly_prod_data$quantity), 2)
p2 <- monthly_prod_data |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "steelblue4") +
geom_vline(xintercept = m, col = "red", size = 1) +
theme_minimal() +
scale_y_continuous(breaks = seq(0, 2500, by = 500), limits = c(0, 2500)) +
# scale_x_binned(breaks = breaks_width(1, 0), limits = c(0,NA)) +
labs(title = paste0("Monthly sales by product ID, median = ",m))Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.print(p2)
The tsibble package has a very useful function to fill in implicit gaps for time series. We could pivot the data wider by product_id, getting a bunch of NA records, and replace them with zero, but this function is convenient and efficient. It uses the first recorded sale as the start, and the last recorded sale as the end, and fills in all gaps with zero
Replacing the implicit gaps with zeros reveals a major challenges with our data, 84% of our records are zeros!
p2b <- data |>
group_by(month, product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
arrange(month) |>
as_tsibble(key = product_id, index = month) |>
fill_gaps(quantity = 0L) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "firebrick4") +
theme_minimal() +
labs(title = "Monthly sales by product ID, implicit gaps removed")
p2b
kable(
data |>
group_by(month, product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
arrange(month) |>
as_tsibble(key = product_id, index = month) |>
fill_gaps(quantity = 0L) |>
as_tibble() |> group_by(quantity) |>
summarize(.count = n()) |>
mutate(
prop = round((.count / sum(.count) * 100), digits = 1)
) |>
top_n(10, prop),
caption = "Frequency of monthly quantity sold by product ID (top 10)"
)| quantity | .count | prop |
|---|---|---|
| 0 | 47190 | 83.5 |
| 1 | 784 | 1.4 |
| 2 | 2062 | 3.7 |
| 3 | 2104 | 3.7 |
| 4 | 1101 | 1.9 |
| 5 | 1172 | 2.1 |
| 6 | 597 | 1.1 |
| 7 | 588 | 1.0 |
| 8 | 281 | 0.5 |
| 9 | 262 | 0.5 |
One method of dealing with intermittent data is temporal aggregation. Instead of looking at just monthly sales, lets look at quarterly and annual.
Data at the quarterly level reduces the number of zeros from over ~47K to almost 12K, and at the annual level we only have around 700.
data |>
group_by(quarter = yearquarter(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = quarter) |>
fill_gaps(quantity = 0L) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "firebrick4") +
theme_minimal() +
labs(title = "Quarterly sales by product ID, implicit gaps removed")
data |>
group_by(year = year(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = year) |>
fill_gaps(quantity = 0L) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "firebrick4") +
theme_minimal() +
labs(title = "Annual sales by product ID, implicit gaps removed")
Maybe some years are behaving differently? For example, as products go obsolete their future values are zero. Unfortunately we don’t have any product management data included so we’ll have to make inferences ourselves.
Breaking it up by year we see that in both 2014 and 2017 we had very few products with zero sales, but in 2015 and 2016 we have many. Why is this?
Our fill gaps function STARTS at the first record, and ENDS on the last record. So for example if a product was sold in only 2014 and 2015, the fill gaps function wont fill in zeros beyond the last sale in 2015.
data |>
group_by(quarter = yearquarter(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = quarter) |>
fill_gaps(quantity = 0L) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "firebrick4") +
theme_minimal() +
labs(title = "Quarterly sales by product ID, implicit gaps removed") +
facet_wrap(vars(year(quarter)))
data |>
group_by(year = year(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = year) |>
fill_gaps(quantity = 0L) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "firebrick4") +
theme_minimal() +
labs(title = "Annual sales by product ID, implicit gaps removed") +
facet_wrap(vars(year))
While our charts above showed very few zeros at the annual level for 2017, when we use another method to just generate a list of products with no sales in 2017, there were actually 337 products. Our data is even MORE intermittent than the charts show!
Ideally we would have information on product life cycle and a process to determine when forecast history starts, when it ends, when replacement products took over, etc. But for this document I’m going to assume all products were available for the full range, and fill in zeros for any period that has no sales.
At this point we have enough information to conclude that forecasting at the product_id level is not going to be pretty. If it was required for the business we could try temporal aggregation, and I come back to address that in section 6.
# kable(
# data |> group_by(product_id, year = year(month)) |>
# summarize(qty = sum(quantity), .groups = "drop") |>
# group_by(year) |>
# summarize(product_count = n(), .groups = "drop"),
# caption = "Unique product count by year"
# )
kable(
nrow(data |> group_by(product_id, year = year(month)) |>
summarize(qty = sum(quantity), .groups = "drop") |>
pivot_wider(names_from = year, values_from = qty) |>
filter(is.na(`2017`))),
caption = "Number of products with no sales in 2017"
)| x |
|---|
| 337 |
# kable(
# data |>
# group_by(product_id) |>
# summarize(total_qty = sum(quantity), .groups = "drop") |>
# top_n(9, total_qty) |>
# arrange(desc(total_qty)),
# caption = "Top 10 products, qty sold 2014-2017"
# )Below are our charts from before with the full range filled in, no reason to drill down to month since will look worse than quarterly, which already looks rough.
data |>
group_by(quarter = yearquarter(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = quarter) |>
fill_gaps(quantity = 0L, .full = TRUE) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "orange2") +
theme_minimal() +
labs(title = "Quarterly sales by product ID, implicit gaps removed full range") +
facet_wrap(vars(year(quarter)))
data |>
group_by(year = year(month), product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = product_id, index = year) |>
fill_gaps(quantity = 0L, .full = TRUE) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, color = "black", fill = "orange2") +
theme_minimal() +
labs(title = "Annual sales by product ID, implicit gaps removed full range") +
facet_wrap(vars(year))
One more view just for fun. This chart takes the first recorded sale of each product and counts their debut by month. It gives us an idea of when products get their first sale, but for course we’re missing all years < 2014 so this picture is incomplete. A very handy function if you have full history, and your company doesn’t track product life cycle information rigorously.
data[match(unique(data$product_id), data$product_id), ] |>
ggplot(aes(x = month)) +
geom_histogram(bins = 48, color = "black", fill = "steelblue4") +
theme_minimal() +
labs(
title = "First recorded sales",
y = "Number of products"
)
Here’s another view (made in hindsight after finishing the project) looking at total sales in 2014-2017, split up by what year they recorded their first sale (start_year on the x axis). What this is getting at is that there is a new product aspect to our data that will cause issues when forecasting. More on this in section 6 when we forecast at the product level.
starts <- data[match(unique(data$product_id), data$product_id), ] |>
select(month, product_id)
data |>
group_by(product_id,year = year(month)) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
left_join(starts, join_by(product_id)) |>
group_by(product_start_year = as.factor(year(month)), year) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(product_start_year, quantity)) +
geom_col(aes(fill = product_start_year)) +
theme_minimal() +
facet_wrap(vars(year), nrow = 1) +
scale_y_continuous(breaks = seq(0,10000, by = 500)) +
labs(title = "Total quantity by year") + theme(axis.text.x = element_text(angle = 90))
3.2.2 Segment
Forecasting for a series of counts (for example units sold) is a complex issue on its own. Forecast values must be integers greater than or equal to zero, and most forecasting models are usually not built to accommodate this. To address the positive integer issue, log transformations can be used when creating a model to force positive forecast values, but as long as the historical volumes are high enough this isn’t an issue.
With that in mind we’re going to look at the different ways to aggregate the data. First being segment as it’s the least familiar to me. The segments aren’t defined and I’m not sure what separates consumer from home office, but I would expect the corporate segment to have a different demand pattern than consumer, and may have different weights at the category/sub_category level.
The first glance at segment sales over time is a little messy to read. It seems like consumer has a different seasonality, but it’s hard to tell because of the different overall sales volume.
data |>
group_by(segment, month) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = segment)) +
geom_line() +
theme_minimal() +
labs(
title = "Total Sales by Month",
x = NULL
) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
for (target_category in unique(data$category)) {
print(data |>
filter(category == target_category) |>
group_by(segment, month) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = segment)) +
geom_line() +
theme_minimal() +
labs(
title = paste(target_category, "Sales by Month"),
x = NULL
)) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
}
The sales by category seem suspiciously similar across the total and all categories, but it’s hard to see for sure at the monthly level. Adding up the total sales history, and viewing the proportions by category shows that while total sales are different by segment, they are split proportionally the same by category! Really unusual and probably a function of the data generation process for this set.
These charts have some interactivity I was experimenting with, so you can hover your mouse and click.
for (target_category in unique(data$category)) {
p14 <- data |>
filter(category == target_category) |>
group_by(segment, month) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = segment)) +
geom_line() +
labs(title = target_category) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") +
labs(x = NULL)
print(plotly::ggplotly(p14))
}
p15 <- data |>
group_by(segment, category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = segment, y = quantity, fill = category)) +
geom_col() +
theme(legend.position = "none") +
theme_minimal() +
theme(legend.position = "none")
p16 <- data |>
group_by(segment, category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = segment, y = quantity, fill = category)) +
geom_col(position = "fill") +
theme_minimal() +
labs(y = "proportion")
plotly::subplot(plotly::ggplotly(p15), plotly::ggplotly(p16))Diving deeper into the sub_category level, we see the same issue.
p16 <- data |>
group_by(segment, sub_category, category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = segment, y = quantity, fill = sub_category)) +
geom_col(position = "fill") +
theme_minimal() +
facet_wrap(vars(category)) +
theme(axis.text.x = element_text(
vjust = 0.5,
angle = 45
)) +
labs(x = NULL, y = "proportion")
plotly::ggplotly(p16)Maybe they have different seasonality and trend components to them? Chart below is a decomposition of these two.
p8 <- data |>
group_by(month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = segment, index = month) |>
model(STL(quantity)) |>
fabletools::components() |>
ggplot(aes(x = month, y = season_year, color = segment)) +
geom_line() +
theme_minimal() +
theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
) +
labs(title = "Trend and seasonality decomposition") +
theme(legend.position = "none")
p9 <- data |>
group_by(month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = segment, index = month) |>
model(STL(quantity)) |>
fabletools::components() |>
ggplot(aes(x = month, y = trend, color = segment)) +
geom_line() +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") +
labs(x = NULL)
print(p8 / p9)
But is this seasonality that we would want to capture through hierarchical forecasting? or is it just noise, hard to tell from the time series chart. A decomposition by category as well will help us see that the segments do actually appear to behave a little differently. While the trend is very similar, the seasonality is clearly different. Consumer has a stronger holiday season in November and December, with a trough in January and February. I’m curious why corporate and consumer share a spike in September, as I thought that for consumer it would be back to school sales but regardless there is a seasonality component here we might want to attempt to capture in our forecast. But as a caveat, there are different total volumes, so it’s difficult to compare by sight.
for (target_category in unique(data$category)) {
p8 <- data |>
filter(category == target_category) |>
group_by(month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = segment, index = month) |>
fill_gaps(quantity = 0) |>
model(STL(quantity)) |>
fabletools::components() |>
ggplot(aes(x = month, y = season_year, color = segment)) +
geom_line() +
labs(title = paste(target_category, "Decomposition")) +
theme_minimal() +
theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
) +
theme(legend.position = "none")
p9 <- data |>
filter(category == target_category) |>
group_by(month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = segment, index = month) |>
fill_gaps(quantity = 0) |>
model(STL(quantity)) |>
fabletools::components() |>
ggplot(aes(x = month, y = trend, color = segment)) +
geom_line() +
theme_minimal() +
theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
p10 <- data |>
filter(category == target_category) |>
group_by(month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = segment, index = month) |>
fill_gaps(quantity = 0) |>
model(STL(quantity)) |>
fabletools::components() |>
ggplot(aes(x = month, y = remainder, color = segment)) +
geom_line() +
theme_minimal() +
tsibble::scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") +
labs(x = NULL) +
theme(legend.position = "none")
print(p8 / p9 / p10)
}
To check more rigorously, we can compare the seasonal strength and trend without being weight by total quantity. There are some differences in trend and seasonality, but whether these differences will give us better forecast accuracy is unclear.
data |>
group_by(category, month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = c("segment", "category"), index = month) |>
fill_gaps(quantity = 0) |>
features(quantity, feat_stl) |>
ggplot(aes(x = trend_strength, y = seasonal_strength_year)) +
geom_point(aes(color = segment), size = 2) +
geom_text(aes(label = category), vjust = "inward", hjust = "inward") +
theme_minimal() +
coord_fixed() +
labs(title = "Seasonal strength vs trend strength") +
xlab(label = "Trend strength") +
ylab(label = "Seasonal strength") +
theme(legend.position = "right")
Not quite as impact when we zoom out however
data |>
group_by(category, month, segment) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = c("segment", "category"), index = month) |>
fill_gaps(quantity = 0) |>
features(quantity, feat_stl) |>
ggplot(aes(x = trend_strength, y = seasonal_strength_year)) +
geom_point(aes(color = segment), size = 2) +
#geom_text(aes(label = category), vjust = "inward", hjust = "inward") +
theme_minimal() +
coord_fixed() +
labs(title = "Seasonal strength vs trend strength") +
xlab(label = "Trend strength") +
ylab(label = "Seasonal strength") +
theme(legend.position = "right") +
scale_x_continuous(limits = c(0, 1)) +
scale_y_continuous(limits = c(0, 1))
3.2.3 Category and Sub Category
Category and sub_category are have a parent/child relationship and products are logically nested based on what type they are. Since we’ve ruled out using forecasting at the product_id level for now, we’ll look into sub_category as our lowest level of aggregation for forecasting, and get our first look at sales over time.
Its clear that we don’t have the same issue of intermittent demand like with product_id. Categories have clear trend and seasonality, and sub_categories have enough volume to generate a forecast.
for (target_category in unique(data$category)) {
p99 <- data |>
filter(category == target_category) |>
group_by(month) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity)) +
geom_line(color = "blue3") +
labs(
title = paste(target_category, "Total"),
x = NULL
) +
theme_minimal() +
tsibble::scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
p14 <- data |>
filter(category == target_category) |>
group_by(month, sub_category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity)) +
geom_line(color = "blue3") +
facet_wrap(vars(sub_category), scales = "free_y") +
labs(
title = paste(target_category, "by Sub Category"),
x = NULL
) +
theme_minimal() +
tsibble::scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
print(p99)
print(p14)
}
51 total unique combinations of category, sub_category, and segment.
data |>
group_by(category, sub_category, segment) |>
summarize(.groups = "drop")# A tibble: 51 × 3
category sub_category segment
<chr> <chr> <chr>
1 Furniture Bookcases Consumer
2 Furniture Bookcases Corporate
3 Furniture Bookcases Home Office
4 Furniture Chairs Consumer
5 Furniture Chairs Corporate
6 Furniture Chairs Home Office
7 Furniture Furnishings Consumer
8 Furniture Furnishings Corporate
9 Furniture Furnishings Home Office
10 Furniture Tables Consumer
# … with 41 more rowsChart below shows the distribution of quantities sold per month using category/sub_category and category/sub_category/segment. At the lowest level (category/sub_category/segment) our data starts getting a lot of zeros, and drilling this far down may lead to less accurate forecasts due to this. But we can prepare models for both with segment and without, to see which performs better.
data |>
group_by(category, sub_category, month) |>
summarize(
quantity = sum(quantity), .groups = "drop"
) |>
as_tsibble(
key = c("category", "sub_category"),
index = month
) |>
fill_gaps(quantity = 0, .full = end()) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, fill = "blue4") +
theme_minimal() +
labs(title = "Category/Sub_category")
data |>
group_by(category, sub_category, segment, month) |>
summarize(
quantity = sum(quantity), .groups = "drop"
) |>
as_tsibble(
key = c("category", "sub_category", "segment"),
index = month
) |>
fill_gaps(quantity = 0, .full = end()) |>
as_tibble() |>
ggplot(aes(x = quantity)) +
geom_histogram(binwidth = 1, fill = "blue4") +
theme_minimal() +
labs(title = "Category/Sub_category/Segment")
4 Develop Models
4.1 Model including segment
Clean data and fill gaps with zero. Code below sets up our history, including segment, as a time series object
range(data$month)<yearmonth[2]>
[1] "2014 Jan" "2017 Dec"fc_data_wsegment <- data |>
group_by(category, sub_category, segment, month) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = c("category", "sub_category", "segment"), index = month)
nrow(has_gaps(fc_data_wsegment) |> filter(.gaps == 1))[1] 40fc_data_wsegment <- fc_data_wsegment |>
fill_gaps(quantity = 0, .full = end())
nrow(has_gaps(fc_data_wsegment) |> filter(.gaps == 1))[1] 0Set up level of aggregation
Category and sub_category have an explicit parent/child relationship. Segment on the other hand will be treated as a group.
fc_data_wsegment <- fc_data_wsegment |>
aggregate_key(((category / sub_category) * segment), quantity = sum(quantity))Subset training data by filtering out 2017
train_data_wsegment <- fc_data_wsegment |>
filter(year(month) < 2017)Fit model to training data
It’s almost unanimous that combining different forecasting models causes less error, so for this approach I use a model combing ETS and ARIMA with equal weights. Rather than modeling and viewing the residuals of each model and fine tuning, we’ll allow the fable package to select the model with the residuals at each level of aggregation. The model selection we’ll do in this document is is whether to use the bottoms up or any min_trace solution.
Reconciliation method selected for this forecast are “mint_shrink” and “ols”, more detail on this at: mint.pdf (robjhyndman.com). In my past experience, mint_shrink has errored out due to low item weight, so we’ll include both just in case.
Ideally we would manually fit a model to each series and adjust trend and seasonality based on input from different departments, this is where enterprise level business forecasting software comes in handy! It can be done in r, but fine tuning 51 different models will make this a little to lengthy.
fit_wsegment <- train_data_wsegment |>
model(
arima = ARIMA(quantity),
ets = ETS(quantity)
) |>
mutate(combo = (arima + ets) / 2) |>
select(-ets, -arima) |>
reconcile(
bu = bottom_up(combo),
shrink = min_trace(combo, method = "mint_shrink"),
ols = min_trace(combo, method = "ols")
)Generate a forecast from our model for the next year (2017) which we can test against our actuals.
fc_wsegment <- fit_wsegment |>
forecast(h = "1 year")Here’s a quick pile of forecast visualizations at different levels compared to our 2017 data. “bu” is the bottom’s up approach, simply summing up the forecasts of the lower levels to get the next level of aggregation “combo” is the direct forecast at a given level of aggregation using the arima and ets models weighted evenly “ols” is the ols method “shrink” is the mint_shrink method. Clearly the shrink and ols are outperforming bottom up, but we’ll check that using MASE further down.
for (target_category in unique(data$category)) {
pfc1 <- fc_wsegment |>
filter(category == target_category,
is_aggregated(segment),
!is_aggregated(sub_category)) |>
autoplot(fc_data_wsegment, alpha = .8, level = NULL) +
facet_wrap(vars(sub_category), scales = "free_y") +
labs(title = target_category) +
theme_minimal() +
theme(axis.text.x = element_text(vjust = 0.5,angle = 45)) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
pfc2 <- fc_wsegment |>
filter(category == target_category,
is_aggregated(segment),
is_aggregated(sub_category)) |>
autoplot(fc_data_wsegment, alpha = .8, level = NULL) +
labs(title = target_category) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
print(pfc2)
print(pfc1)
}
fc_wsegment |>
filter(is_aggregated(segment), is_aggregated(sub_category), is_aggregated(category)) |>
autoplot(fc_data_wsegment, alpha = .8, level = NULL) +
facet_wrap(vars(.model)) +
theme_minimal() +
labs(title = "Total") +
theme(axis.text.x = element_text(
vjust = 0.5,
angle = 45
)) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
fc_wsegment |>
filter(!is_aggregated(segment), is_aggregated(sub_category), is_aggregated(category)) |>
autoplot(fc_data_wsegment, alpha = .8, level = NULL) +
facet_wrap(vars(segment)) +
theme_minimal() +
labs(title = "Segment") +
theme(axis.text.x = element_text(
vjust = 0.5,
angle = 45
)) +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
4.2 Model without segment
And here we do all the same without including segment, just category/sub_category.
Fill gaps without segment.
fc_data_wosegment <- data |>
group_by(category, sub_category, month) |>
summarize(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(key = c("category", "sub_category"), index = month)
nrow(has_gaps(fc_data_wosegment) |> filter(.gaps == 1))[1] 6fc_data_wosegment <- fc_data_wosegment |>
fill_gaps(quantity = 0, .full = end())
nrow(has_gaps(fc_data_wosegment) |> filter(.gaps == 1))[1] 0Set up aggregation.
fc_data_wosegment <- fc_data_wosegment |> aggregate_key(category / sub_category, quantity = sum(quantity))Subset training data, filtering out 2017.
train_data_wosegment <- fc_data_wosegment |> filter(year(month) < 2017)Fit model to training data.
fit_wosegment <- train_data_wosegment |>
model(
arima = ARIMA(quantity),
ets = ETS(quantity)
) |>
mutate(combo = (arima + ets) / 2) |>
select(-ets, -arima) |>
reconcile(
bu = bottom_up(combo),
shrink = min_trace(combo, method = "mint_shrink"),
ols = min_trace(combo, method = "ols")
)Generate a forecast from our model.
fc_wosegment <- fit_wosegment |>
forecast(h = "1 year")Forecast visualization.
for (target_category in unique(data$category)) {
pfc1 <- fc_wosegment |>
filter(category == target_category, !is_aggregated(sub_category)) |>
autoplot(fc_data_wosegment, alpha = .8, level = NULL) +
facet_wrap(vars(sub_category), scales = "free_y") +
labs(title = target_category) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
pfc2 <- fc_wosegment |>
filter(category == target_category, is_aggregated(sub_category)) |>
autoplot(fc_data_wosegment, alpha = .8, level = NULL) +
labs(title = target_category) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y")
print(pfc2)
print(pfc1)
}
fc_wosegment |>
filter(is_aggregated(sub_category), is_aggregated(category)) |>
autoplot(fc_data_wosegment, alpha = .8, level = NULL) +
facet_wrap(vars(.model)) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") +
labs(title = "Total")
5 Forecast Accuracy comparison
As expected, the bottom up approach gets less accurate when we add segmentation because each series has less information
Directly forecasting the total was the most accurate at that level
min trace solutions created forecasts with the lowest MASE for category and sub_category. Minor differences between ols and mint shrink
Segmentation in general only had a minor impact to MASE
The decisions to make here are which reconciliation method to use, and whether to include segmentation in our model. While segmentation didn’t impact forecast performance too much, I’m going to choose to create the final forecast using it. There are some theoretical benefits to the business for having forecasts based on region and it isn’t making our test forecast worse, so we’ll assume there is a hypothetical business need for it. I chose the min_trace reconciliation method using mint_shrink because it created slightly more accurate forecasts at the lowest level.
fa1 <- fc_wsegment |>
filter(is_aggregated(category), is_aggregated(sub_category), is_aggregated(segment)) |>
accuracy(
data = fc_data_wsegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "yes", level = "total")
fa2 <- fc_wosegment |>
filter(is_aggregated(category), is_aggregated(sub_category)) |>
accuracy(
data = fc_data_wosegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "no", level = "total")
fa3 <- fc_wsegment |>
filter(!is_aggregated(category), is_aggregated(sub_category), is_aggregated(segment)) |>
accuracy(
data = fc_data_wsegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "yes", level = "category")
fa4 <- fc_wosegment |>
filter(!is_aggregated(category), is_aggregated(sub_category)) |>
accuracy(
data = fc_data_wosegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "no", level = "category")
fa5 <- fc_wsegment |>
filter(!is_aggregated(category), !is_aggregated(sub_category), is_aggregated(segment)) |>
accuracy(
data = fc_data_wsegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "yes", level = "sub_category")
fa6 <- fc_wosegment |>
filter(!is_aggregated(category), !is_aggregated(sub_category)) |>
accuracy(
data = fc_data_wosegment,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(include_segment = "no", level = "sub_category")
accuracytable <- bind_rows(fa1, fa2, fa3, fa4, fa5, fa6)
accuracytable <- accuracytable |> pivot_wider(names_from = .model, values_from = mase)
kable(accuracytable,
caption = "MASE comparison \n
bu = Bottom-up, \n
combo = base forecast at that level of aggregation \n
ols = min trace ols reconciliation method \n
shrink = min trace shrink reconciliation method",
digits = 2
)| include_segment | level | bu | combo | ols | shrink |
|---|---|---|---|---|---|
| yes | total | 2.08 | 1.30 | 1.39 | 1.40 |
| no | total | 1.60 | 1.30 | 1.33 | 1.39 |
| yes | category | 1.96 | 1.48 | 1.42 | 1.41 |
| no | category | 1.56 | 1.48 | 1.41 | 1.44 |
| yes | sub_category | 1.30 | 1.22 | 1.16 | 1.14 |
| no | sub_category | 1.22 | 1.22 | 1.18 | 1.17 |
5.1 Final Forecast
Usually I try to avoid “final” because it leads to files names like final_rev2, and final_final, but I think it’s okay for this document. <later edit, I was wrong>
fit_wsegment <- fc_data_wsegment |>
model(
arima = ARIMA(quantity),
ets = ETS(quantity)
) |>
mutate(combo = (arima + ets) / 2) |>
select(-ets, -arima) |>
reconcile(
mint = min_trace(combo, method = "mint_shrink")
)Forecast out one year
fc <- fit_wsegment |>
forecast(h = "1 year")Final forecasts
fc |>
filter(is_aggregated(sub_category), is_aggregated(segment), .model == "mint") |>
autoplot(fc_data_wsegment, level = 90) +
facet_wrap(vars(category), scales = "free_y") +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") +
labs(title = "Category")
for (target_category in unique(data$category)) {
p <- fc |>
filter(
category == target_category, !is_aggregated(sub_category),
is_aggregated(segment),
.model == "mint"
) |>
autoplot(fc_data_wsegment, level = 90) +
facet_wrap(vars(sub_category), scales = "free_y") +
labs(title = target_category) +
theme_minimal() +
scale_x_yearmonth(date_breaks = "1 year", date_labels = "%Y") + theme(axis.text.x = element_text(
vjust = 0.5,
angle = 45
))
print(p)
}
5.2 Clean up the final forecast for export, remove zeros, change to integer values
Example of what exported data would look like, issues rounding leads to incoherent forecasts and removing zeros does as well. Be careful when interpreting and publishing the different levels of aggregation.
final_fcst <- fc |>
filter(.model == "mint") |>
as_tibble() |>
select(-.model) |>
mutate(
date = make_date(year = year(month), month = month(month), day = 1),
`80%` = hilo(quantity, 80),
`90%` = hilo(quantity, 90)
) |>
select(-quantity, -month) |>
unpack_hilo(c(`80%`, `90%`)) |>
relocate(date)
final_fcst <- final_fcst |>
mutate(across(
.cols = .mean:`90%_upper`,
.fns = ~ pmax(., 0)
)) |>
mutate(across(
.cols = .mean:`90%_upper`,
.fns = ~ as.integer(.x)
))
write_csv(final_fcst, file = "final_fcst")
kable(head(final_fcst, n = 10),
caption = "Example of data export \n"
)| date | category | segment | sub_category | .mean | 80%_lower | 80%_upper | 90%_lower | 90%_upper |
|---|---|---|---|---|---|---|---|---|
| 2018-01-01 | Furniture | Consumer | Bookcases | 10 | 1 | 19 | 0 | 22 |
| 2018-02-01 | Furniture | Consumer | Bookcases | 11 | 2 | 20 | 0 | 22 |
| 2018-03-01 | Furniture | Consumer | Bookcases | 15 | 6 | 24 | 3 | 27 |
| 2018-04-01 | Furniture | Consumer | Bookcases | 17 | 7 | 26 | 5 | 29 |
| 2018-05-01 | Furniture | Consumer | Bookcases | 17 | 7 | 26 | 4 | 29 |
| 2018-06-01 | Furniture | Consumer | Bookcases | 19 | 9 | 28 | 6 | 31 |
| 2018-07-01 | Furniture | Consumer | Bookcases | 15 | 6 | 25 | 3 | 28 |
| 2018-08-01 | Furniture | Consumer | Bookcases | 14 | 4 | 24 | 2 | 26 |
| 2018-09-01 | Furniture | Consumer | Bookcases | 22 | 12 | 32 | 9 | 35 |
| 2018-10-01 | Furniture | Consumer | Bookcases | 18 | 9 | 28 | 6 | 31 |
6 Coming back to this project to forecast at product level
I was planning on leaving this analysis to the sub_category level, but I ended up coming back later wanting to create product level forecasts. Turns out I shouldn’t have named that table “final_fcst” after all!
As much as I would love to only forecast high volume series, business forecasting doesn’t usually have that luxury. Creating category and sub_category level forecast can be extremely helpful to many parts of the a business, from trend analysis to allocating sales resources to factory and staff planning, but forecasts are often used for supply chain purchasing purposes which require forecasts at the lowest level of aggregation.
In this section I’m going to do the following and address the issues caused by each
Analyse item level at daily, weekly, monthly, quarterly, and annual to see if we can see trends at any of those levels
- We built a histogram of historical sales quantities at these levels, but we’ll dive deeper in this section
Apply Cronston model which is built for intermittent demand
Disaggregate the Category/Sub_Category forecast, weighted by cronston forecast (come back to this if we have poor performance with reconciliation method) This is a less fancy solution but would work
Apply temporal hierarchies combined with cronston
- We’ll quickly run into issues using the cronston model that conflicts with a “one number” plan/forecast. We’ll see if we can take the best of both worlds and create a model that creates sub_category, category, and total forecasts that make sense, while creating also creating product level forecasts
Convert forecast values to integers
- We can’t just round like before because many if not most forecast values will be less than one at the item level
6.1 Reload the data and clean
We used the “data” variable in some transformations earlier so we can’t use it again. Lets reload it as data2 and clean it up with the same code as before.
data2 <- readxl::read_xls("C:/Users/andre/Documents/r4ds/SampleSuperStore/Sample - Superstore.xls")
# clean up column names, using snake_case as default
data2 <- data2 |> janitor::clean_names()
data2 <- data2 |>
mutate(order_date = as.Date(order_date), ship_date = as.Date(ship_date), .keep = "all") |>
select(-row_id)
data2 <- data2 |>
group_by(
order_date,
product_id,
category,
sub_category
) |>
summarize(quantity = sum(quantity), .groups = "drop")Set up a list 20 of random products to take a look at.
set.seed(99)
product_list <- unique(data2$product_id)
product_sample <- sample(product_list, 10)
product_sample [1] "FUR-BO-10000780" "FUR-FU-10001037" "OFF-PA-10003205" "OFF-PA-10000726"
[5] "FUR-CH-10002965" "TEC-AC-10004568" "OFF-AR-10002956" "FUR-CH-10003774"
[9] "FUR-FU-10003347" "TEC-CO-10001766"Set up our time series table with category, sub_category, and product_id. I chose to drop segment in this case, we already know that our product_id level data is sparse, so no reason to divide it into three groups by segment. Here we take our data; summarize on order_date, category, sub_category, and product_id; coerce to a tsibble, and fill in the gaps for the full range in one step.
data2_ts <- data2 |>
group_by(order_date, category, sub_category, product_id) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = order_date
) |>
fill_gaps(quantity = 0, .full = TRUE)Quick look at the daily data, yep its a mess.
data2_ts |>
filter(product_id %in% product_sample) |>
ggplot(aes(x = order_date, y = quantity)) +
geom_line(aes(color = product_id)) +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "Random 10 product_id",
x = "day"
)
6.2 Check trend/seasonality at different temporal aggregations
Set up weekly, monthly, quarterly, and annual time series, maybe we can get some value from temporal aggregation (aggregating time to these different levels)
data2_ts_weekly <- data2_ts |>
as_tibble() |>
group_by(
week = yearweek(order_date),
product_id,
category,
sub_category
) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = week
)
data2_ts_monthly <- data2_ts |>
as_tibble() |>
group_by(
month = yearmonth(order_date),
product_id,
category,
sub_category
) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = month
)
data2_ts_quarterly <- data2_ts |>
as_tibble() |>
group_by(
quarter = yearquarter(order_date),
product_id,
category,
sub_category
) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = quarter
)
data2_ts_annual <- data2_ts |>
as_tibble() |>
group_by(
year = year(order_date),
product_id,
category,
sub_category
) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = year
)Plot that same random sample… and we see that it’s still essentially random noise, we’re not getting any trend or seasonality. If anything there might
p1 <- data2_ts_weekly |>
filter(product_id %in% product_sample) |>
ggplot(aes(x = week, y = quantity, color = product_id, group = product_id)) +
geom_line() +
theme_minimal() +
theme(legend.position = "none")
p2 <- data2_ts_monthly |>
filter(product_id %in% product_sample) |>
ggplot(aes(x = month, y = quantity, color = product_id)) +
geom_line() +
theme_minimal() +
theme(legend.position = "none")
p3 <- data2_ts_quarterly |>
filter(product_id %in% product_sample) |>
select(product_id, quarter, quantity) |>
ggplot(aes(x = quarter, y = quantity, color = product_id, group = product_id)) +
geom_line() +
theme_minimal() +
theme(legend.position = "none")
p4 <- data2_ts_annual |>
filter(product_id %in% product_sample) |>
select(product_id, year, quantity) |>
ggplot(aes(x = year, y = quantity, color = product_id)) +
geom_line() +
theme_minimal() +
theme(legend.position = "none")
((p1 | p2) / (p3 | p4)) + plot_annotation(title = "random 20 product_id")
How about the top 10 selling products? Still not great, we would expect very clear trend and seasonality from the highest volume, but we’re not getting it. At this point we can assume all product data is intermittent and we can use the Cronston method to forecast them, and forget the idea of using temporal aggregation. We might have some trend at the annual level, but would be relying on four periods of history to generate a forecast.
top10 <- data2_ts |>
as_tibble() |>
group_by(product_id) |>
summarise(quantity = sum(quantity)) |>
arrange(desc(quantity)) |>
top_n(10, quantity) |>
select(product_id)
p1 <- data2_ts_weekly |>
filter(product_id %in% top10$product_id) |>
ggplot(aes(x = week, y = quantity, color = product_id)) +
geom_line() +
theme_minimal() +
theme(legend.position = "none") +
labs(title = "weekly")
p2 <- data2_ts_monthly |>
filter(product_id %in% top10$product_id) |>
ggplot(aes(x = month, y = quantity, color = product_id)) +
geom_line() +
theme(legend.position = "none") +
theme_minimal() +
labs(title = "monthly") +
theme(legend.position = "none")
p3 <- data2_ts_quarterly |>
filter(product_id %in% top10$product_id) |>
ggplot(aes(x = quarter, y = quantity, color = product_id)) +
geom_line() +
theme(legend.position = "none") +
theme_minimal() +
labs(title = "quarterly") +
theme(legend.position = "none")
p4 <- data2_ts_annual |>
filter(product_id %in% top10$product_id) |>
ggplot(aes(x = year, y = quantity, color = product_id)) +
geom_line() +
theme(legend.position = "none") +
theme_minimal() +
labs(title = "annual") +
theme(legend.position = "none")
(p1 | p2) / (p3 | p4) + plot_annotation(title = "Top 10 product_id")
6.3 Forecast with Crontson for intermittent demand
I defaulted here to monthly data, which from my experience is the typical period for business forecasting. Then fit a model using the Cronston method, and forecast 12 months, using 2014-2016 to fit the model, and 2017 to get forecast accuracy metrics. The test model generates 297 null models, about 16% of our products are erroring out! This really isn’t a problem now, since the null models will generate a NA forecast, but it will be a problem if we chose to reconcile it using the reconcile function in the fable package. The null models are likely because the history for that product only has one period of sales, and we need at least two periods to create a model.
monthly_cronston_fit_test <- data2_ts_monthly |>
filter(year(month) < 2017) |>
model(.model = CROSTON(quantity))Warning: 297 errors (1 unique) encountered for .model
[297] At least two non-zero values are required to use Croston's method.monthly_cronston_fc_test <- monthly_cronston_fit_test |>
forecast(h = 12)We can test this assumption by filtering to the null models from our model table, and then filter our sales data by that product list with sales > 0. Length of this table is 297, meaning that this is indeed the cause of our errors. Null models will create NA forecast values, which isn’t a problem for now, we can just replace those with zeros, and in a business process they would be marked for review.
We’ll create the list of products and an object we can filter for later
monthly_cronston_nulls_test <- monthly_cronston_fit_test |>
mutate(null = is_null_model(.model)) |>
filter(null == 1)
length(unique(monthly_cronston_nulls_test$product_id))[1] 297Lets take a look at total quantities over time, and yikes! These products have a lot more sales in 2017. This will make things goofy comparing accuracy to our previous method (using only category and sub_category). Our category/sub_category forecast was under forecasting, good reason to think these products may be why we had so much growth in 2017, and would be a good focus for further analysis (not here though)
data2_ts_monthly |> filter(product_id %in% monthly_cronston_nulls_test$product_id) |>
index_by(month) |>
summarise(quantity = sum(quantity)) |>
ggplot(aes(x = month, y = quantity)) +
geom_line() +
theme_minimal()
When we add up the cronston forecasts, we quickly see a problem, they’re flat! Cronston method produces forecasts with no trend, no seasonality, and no confidence range. Not very helpful if using a bottom up approach!
asdf <- data2_ts_monthly |>
as_tibble() |>
mutate(type = "history") |>
select(product_id, category, sub_category, type, month, quantity) |>
group_by(category, month, type) |>
summarize(quantity = sum(quantity), .groups = "drop")
asdf2 <- monthly_cronston_fc_test |>
filter(!(product_id %in% monthly_cronston_nulls_test$product_id)) |>
as_tibble() |>
mutate(type = "forecast") |>
group_by(category, month, type) |>
summarize(quantity = sum(.mean), .groups = "drop")
asdf3 <- rbind(asdf,asdf2)
for(target_category in unique(asdf3$category)) {
p1 <- asdf3 |> filter(category == target_category) |>
ggplot(aes(month, quantity, col = type)) +
geom_line() +
theme_minimal() +
labs(title = target_category) +
scale_color_manual(values = c("red","black"))
print(p1)
}
6.3.1 Cronston and proportional split
This section is currently under development. Idea is to do a napkin math method of proportionally splitting item forecasts based on category/sub_category forecast. This would be a great option to create a forecast without worrying about the rigid structure of this forecasting package
https://robjhyndman.com/papers/foresight.pdf
Lets generate MASE metrics for our 2017 forecast, taking the mean by category
As Marcel the Shell famously said, “compared to what??”; the metrics aren’t very meaningful when there’s nothing to compare to yet. But we can tell that accuracy is not really great, forecasted values are worse than the one step ahead naive forecast. but this makes sense, considering that the most accurate forecast would probably be zeros across the board.
# monthly_cronston_fc_accuracy_test <-
# monthly_cronston_fc_test |>
# filter(!(product_id %in% monthly_cronston_nulls_test$product_id)) |>
# accuracy(
# data = data2_ts_monthly |>
# filter(!(product_id %in% monthly_cronston_nulls_test$product_id)
# ),
# measures = list(mase = MASE)
# )
#
# kable(monthly_cronston_fc_accuracy_test |>
# group_by(category) |>
# summarise(MASE = mean(mase)), caption = "cronston test accuracy")Also, compare the 2017 forecasts(not our test) to our previous test forecast by category/sub_category. The total by category actually roll up fairly close, but we quickly see a the problem with the Cronston method when we chart by month. The forecast is flat when we roll it up! This is always the case when using this model, its effective for intermittent demand, but creates a flat forecast (no trend or seasonality). This item level forecast might be helpful for supply planning on intermittent products, but if we roll it up it doesn’t make much sense ***** come back to, need to discuss models that dropped out that make comparison difficult, maybe just view v actuals*****
# monthly_cronston_fc_test |>
# as_tibble() |>
# group_by(category) |>
# summarise(quantity = sum(.mean, na.rm = T))
#
# fc_wosegment |>
# as_tibble() |>
# filter(!is_aggregated(category),
# !(as.character(product_id) %in% monthly_cronston_nulls_test$product_id),
# .model == "ols"
# ) |>
# group_by(category) |>
# summarise(quantity = sum(.mean, na.rm = T))
#
# fc_method_compare <- maybe_works_fit2_rec_fc_test |>
# filter(
# !is_aggregated(product_id),
# !(as.character(product_id) %in% monthly_cronston_nulls_test$product_id),
# .model == "mint"
# ) |>
# as_tibble() |>
# group_by(month, category) |>
# summarise(aggregated_fc = sum(.mean)) |>
# left_join(monthly_cronston_fc_test |> as_tibble() |> group_by(month, category) |> summarise(cronston_fc = sum(.mean, na.rm = T)), join_by("month" == "month", "category" == "category")) |>
# left_join(
# data2_ts_monthly |> filter(
# !(product_id %in% monthly_cronston_nulls_test$product_id),
# year(month) == 2017) |> as_tibble() |>
# group_by(month, category) |>
# summarise(actuals = sum(quantity)),
# join_by("month" == "month", "category" == "category")
# )
#
#
#
#
# fc_method_compare |>
# pivot_longer(cols = c(aggregated_fc, cronston_fc, actuals), names_to = "method", values_to = "quantity") |>
# mutate(category = as.character(category)) |>
# ggplot(aes(x = month, y = quantity)) +
# geom_line(aes(color = method)) +
# facet_wrap(vars(category), scales = "free", ncol = 1) +
# theme_minimal() +
# labs(title = "2017 sales")One way to fix this is to take the proportions of our product_id level forecast, and apply those across sub_categories using our previous reconciled forecast at that level ****come back to this too****.
# cronstonfc_prop_test <- monthly_cronston_fc_test |>
# as_tibble() |>
# filter(!(product_id %in% monthly_cronston_nulls_test))
# select(product_id, category, sub_category, .mean) |>
# mutate(.mean = replace_na(.mean, 0)) |>
# group_by(product_id, category, sub_category) |>
# summarise(qty = sum(.mean), .groups = "drop") |>
# group_by(sub_category) |>
# mutate(prop = qty / sum(qty))
# data2_ts |> as_tibble() |> group_by(category, sub_category, product_id) |>
# summarise(quantity = sum(quantity), .groups = "drop") |>
# group_by(category) |> mutate(cat_prop = quantity/sum(quantity)) |> ungroup() |>
# group_by(sub_category) |> mutate(sub_cat_prop = quantity/sum(quantity)) |> ungroup() |>
# mutate(total_prop = quantity/sum(quantity)) |> select(-quantity)Check to make sure totals are rolling up correctly when using this method, compared to our prior forecast ****come back to this too****.
# sum up final forecast and drop segment for comparison, drop range because it wont be relevant at the product level
# final_fc_for_props <- final_fcst |>
# filter(is_aggregated(segment)) |>
# select(-segment) |>
# select(date:.mean) |>
# mutate(category = as.character(category), sub_category = as.character(sub_category))
#
# asdf <- final_fc_for_props |>
# right_join(cronstonfc_prop, join_by("category" == "category", "sub_category" == "sub_category")) |>
# mutate(item_fc = (.mean * prop))
#
# asdf |>
# group_by(category) |>
# summarise(qty = sum(item_fc))
# final_fcst |>
# filter(!is_aggregated(segment), !is_aggregated(sub_category)) |>
# group_by(category) |>
# summarise(qty = sum(.mean))6.4 Can we use the hierarchical forecast approach?
Unfortunately the fable package doesn’t yet have a ‘case when’ option for setting models, so while it’s a little time consuming, we’ll set an ETS model and cronston model to all products and levels of aggregation. remember those 97 products that errored out before? Now they’ll become an issue, because the reconciling approach cannot accept null models. Lets drop those out before fitting the model.
Lets set up our data to the monthly level with the item hierarchy (category/sub_category/product_id) and go through this approach, dropping the 297 products that only have one observed sale period.
data2_ts_monthly_ag_test <- data2_ts_monthly |>
filter(year(month) < 2017,
!(product_id %in% monthly_cronston_nulls_test$product_id)) |>
aggregate_key(category / sub_category / product_id, quantity = sum(quantity))6.4.1 Set up model
Fit arima/ets combo and cronston forecast.
ttl_ag_fit_test <-
data2_ts_monthly_ag_test |> model(
arima = ARIMA(quantity),
ets = ETS(quantity),
cronston = CROSTON(quantity)
)
ttl_ag_fit_test <- ttl_ag_fit_test |>
mutate(combo = (arima + ets) / 2) |>
select(-ets, -arima)The mable (model table) created is a complicated table bound by the aggregation hierarchy and date index. Each record in the table is actually hundreds of nested lists. It look a lot of head scratching and forum searching, but I found that if we split them by level of aggregation and use bind_rows, we can use the cronston forecast at the product level, and our arima/ets combo for the sub_category and category level models.
ags <- ttl_ag_fit_test |> filter(is_aggregated(product_id)) |> select(-cronston) |> rename(.model = combo)
prods <- ttl_ag_fit_test |> filter(!is_aggregated(product_id)) |> select(-combo) |> rename(.model = cronston)
ags# A mable: 21 x 4
# Key: category, sub_category, product_id [21]
category sub_category product_id .model
<chr*> <chr*> <chr*> <model>
1 Furniture Bookcases <aggregated> <COMBINATION>
2 Furniture Chairs <aggregated> <COMBINATION>
3 Furniture Furnishings <aggregated> <COMBINATION>
4 Furniture Tables <aggregated> <COMBINATION>
5 Furniture <aggregated> <aggregated> <COMBINATION>
6 Office Supplies Appliances <aggregated> <COMBINATION>
7 Office Supplies Art <aggregated> <COMBINATION>
8 Office Supplies Binders <aggregated> <COMBINATION>
9 Office Supplies Envelopes <aggregated> <COMBINATION>
10 Office Supplies Fasteners <aggregated> <COMBINATION>
# … with 11 more rowsprods# A mable: 1,565 x 4
# Key: category, sub_category, product_id [1,565]
category sub_category product_id .model
<chr*> <chr*> <chr*> <model>
1 Furniture Bookcases FUR-BO-10000330 <croston>
2 Furniture Bookcases FUR-BO-10000362 <croston>
3 Furniture Bookcases FUR-BO-10000468 <croston>
4 Furniture Bookcases FUR-BO-10000711 <croston>
5 Furniture Bookcases FUR-BO-10000780 <croston>
6 Furniture Bookcases FUR-BO-10001337 <croston>
7 Furniture Bookcases FUR-BO-10001519 <croston>
8 Furniture Bookcases FUR-BO-10001608 <croston>
9 Furniture Bookcases FUR-BO-10001798 <croston>
10 Furniture Bookcases FUR-BO-10001811 <croston>
# … with 1,555 more rowsttl_ag_fit_test_bind <- bind_rows(ags,prods)Set up reconciliation method using ols. I tried using mint_shrink before, but there is a line in the function causing it to error out if it finds an item with too small of a weight.
ttl_ag_fit_test_bind_rec <- ttl_ag_fit_test_bind |>
reconcile(ols = min_trace(.model, method = "ols"))Forecast out 12 months.
ttl_ag_fit_test_fc <- ttl_ag_fit_test_bind_rec |>
forecast(h = 12)Lets visualize how this forecast looks compared to actuals, and it looks pretty reasonable! It looks even more accurate (esp for late 2017) than our sub_category/category forecast, but that’s because we dropped those 297 products that had huge growth in 2017.
data2_ts_monthly_ag_test_full <- data2_ts_monthly |>
filter(!(product_id %in% monthly_cronston_nulls_test$product_id)) |>
aggregate_key(category / sub_category / product_id, quantity = sum(quantity))
ttl_ag_fit_test_fc |> filter(is_aggregated(category)) |>
autoplot(data2_ts_monthly_ag_test_full, level = NULL) +
theme_minimal()
ttl_ag_fit_test_fc |> filter(is_aggregated(sub_category), !is_aggregated(category)) |>
autoplot(data2_ts_monthly_ag_test_full, level = NULL) +
theme_minimal() +
facet_wrap(vars(category), ncol = 3, scales = "free")
ttl_ag_fit_test_fc |> filter(category == "Furniture", !is_aggregated(sub_category), is_aggregated(product_id)) |>
autoplot(data2_ts_monthly_ag_test_full, level = NULL) +
theme_minimal() +
facet_wrap(vars(sub_category), scales = "free")
6.4.2 Compare to 2017 actuals
Lets check the MASE for this new forecast method…and we see that this method is actually hurting our accuracy across all levels, compared to the direct forecast at that level of aggregation. But the benefit is that this product level forecast rolls up to a total forecast that actually makes sense and fits history! Depending on the business needs this is a worthwhile trade off because the accuracy losses are very minimal.
A forecaster should of course have a process to review products with the highest MASE score to determine root cause, but that’s outside my scope for this document.
Another thing this tells us is how impactful those 297 dropped products where to forecast accuracy. This gains in MASE vs our forecast in section 5 aren’t from the detailed product_id level forecasts, but because products with their second sale in 2017 accounted for a significant portion of the growth
-We could rerun our category/sub_category forecast, dropping out those 297 products to get a fair comparison between the methods but this report already takes an hour to render with all the modeling so I’m not going to do it here
fa1 <- ttl_ag_fit_test_fc |>
filter(is_aggregated(category), is_aggregated(sub_category), is_aggregated(product_id)) |>
accuracy(
data = data2_ts_monthly_ag_test_full,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(level = "total")
fa2 <- ttl_ag_fit_test_fc |>
filter(!is_aggregated(category), is_aggregated(sub_category), is_aggregated(product_id)) |>
accuracy(
data = data2_ts_monthly_ag_test_full,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(level = "category")
fa3 <- ttl_ag_fit_test_fc |>
filter(!is_aggregated(category), !is_aggregated(sub_category), is_aggregated(product_id)) |>
accuracy(
data = data2_ts_monthly_ag_test_full,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(level = "sub_category")
fa4 <- ttl_ag_fit_test_fc |>
filter(!is_aggregated(category), !is_aggregated(sub_category), !is_aggregated(product_id)) |>
accuracy(
data = data2_ts_monthly_ag_test_full,
measures = list(mase = MASE)
) |>
group_by(.model) |>
summarize(mase = mean(mase)) |>
mutate(level = "product_id")
accuracytable2 <- bind_rows(fa1, fa2, fa3, fa4)
accuracytable2 <- accuracytable2 |> pivot_wider(names_from = .model, values_from = mase)
kable(accuracytable2,
digits = 2)| level | .model | ols |
|---|---|---|
| total | 0.90 | 0.91 |
| category | 0.98 | 1.01 |
| sub_category | 1.01 | 1.01 |
| product_id | 1.56 | 1.63 |
6.4.3 Forecast 2018
Lets save some computing time and calculate check which products only have one sales record. Generates list of 97, we’ll drop these out at the beginning because we’ll get null models.
data2_ts_monthly |> as_tibble() |> mutate(helper = ifelse(quantity > 0, 1, 0)) |> group_by(product_id) |>
summarize(helper = sum(helper)) |> filter(helper < 2) |> nrow() |> kable(caption = "Number of products with < 2 periods of sales")| x |
|---|
| 97 |
drops <- data2_ts_monthly |> as_tibble() |> mutate(helper = ifelse(quantity > 0, 1, 0)) |> group_by(product_id) |>
summarize(helper = sum(helper)) |> filter(helper < 2) |> select(product_id)Create table dropping out those products.
data2_ts_monthly_wdrops <- data2_ts_monthly |> filter(!product_id %in% drops$product_id)create time series object without those products at monthly level.
data2_ts_monthly_wdrops <- data2_ts |>
as_tibble() |>
filter(!(product_id %in% drops$product_id)) |>
group_by(
month = yearmonth(order_date),
product_id,
category,
sub_category
) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
as_tsibble(
key = c("product_id", "category", "sub_category"),
index = month
)create aggregation.
data2_ts_monthly_wdrops_ag <- data2_ts_monthly_wdrops |>
aggregate_key(category / sub_category / product_id, quantity = sum(quantity))fit ets/arima combo and cronston model. combo will be used for cat thru sub_cat, cronston for product level.
ttl_ag_wdrops_fit <-
data2_ts_monthly_wdrops_ag |> model(
ets = ETS(quantity),
arimia = ARIMA(quantity),
cronston = CROSTON(quantity)
)Warning in sqrt(diag(best$var.coef)): NaNs producedttl_ag_wdrops_fit <- ttl_ag_wdrops_fit |>
mutate(combo = (arimia + ets) / 2) |>
select(-ets, -arimia)splice our mable (model table) by product hierarchy, select appropriate model for each level, rebind it together. dplyer makes it easy to bind matrices together using bind_rows. In the code below ags are the cat thru sub_cat models, prods are the product_id level models.
ags <- ttl_ag_wdrops_fit |>
filter(is_aggregated(product_id)) |>
select(-cronston) |>
rename(.model = combo)
prods <- ttl_ag_wdrops_fit |>
filter(!is_aggregated(product_id)) |>
select(-combo) |>
rename(.model = cronston)
ttl_ag_wdrops_fit <- bind_rows(ags, prods)take our new mable and add reconcile method ols.
ttl_ag_wdrops_fit_rec <-
ttl_ag_wdrops_fit |>
reconcile(ols = min_trace(.model, method = "ols"))forecast out 12 months.
ttl_ag_wdrops_fit_rec_fc <- ttl_ag_wdrops_fit_rec |>
forecast(h = 12)visuals comparing section 5 and section 6 forecasts for 2018
for (target_category in unique(data2$category)) {
p1 <- ttl_ag_wdrops_fit_rec_fc |>
filter(category == target_category,
is_aggregated(sub_category),
is_aggregated(product_id),
.model == "ols") |>
autoplot(data2_ts_monthly_wdrops_ag, level = NULL) +
theme_minimal() +
theme(legend.position = "none") +
labs(title = paste0(target_category," ", "Category/Sub_category/Product_id", " ", "Forecast"))
p2 <- fc |>
filter(category == target_category,
is_aggregated(sub_category), is_aggregated(segment), .model == "mint") |>
autoplot(fc_data_wsegment, level = NULL) +
theme_minimal() +
labs(title = paste0(target_category," ", "Category/Sub_category", " ", "Forecast"))
print(p1/p2)
}
Chug out some more visuals comparing our forecasts for 2018
catsubcat_tibble <- as_tibble(fc) |> filter(!is_aggregated(category), !is_aggregated(segment), !is_aggregated(sub_category),
.model == "mint") |>
mutate(category = as.character(category),
sub_category = as.character(sub_category),
.model = "cat/sub_cat Forecast",
month = yearmonth(month),
quantity = .mean) |>
select(
category, sub_category, .model, month, quantity) |>
group_by(category, sub_category, month, .model) |>
summarise(quantity = sum(quantity), .groups = "drop")
catsubcatprod_tibble <- as_tibble(ttl_ag_wdrops_fit_rec_fc) |>
filter(!is_aggregated(category), !is_aggregated(sub_category), !is_aggregated(product_id),
.model == "ols") |>
mutate(category = as.character(category),
sub_category = as.character(sub_category),
.model = "cat/sub_cat/product_id Forecast",
month = yearmonth(month),
quantity = .mean) |>
select(
category, sub_category, .model, month, quantity) |>
group_by(category, sub_category, month, .model) |>
summarise(quantity = sum(quantity), .groups = "drop")
history <- data2_ts_monthly |>
as_tibble() |>
group_by(month, category, sub_category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
mutate(.model = "history")
catsubcat_tibble <- catsubcat_tibble |> select(month, category, sub_category, quantity, .model)
catsubcatprod_tibble <- catsubcatprod_tibble |> select(month, category, sub_category, quantity, .model)
asdf <- rbind(catsubcat_tibble, catsubcatprod_tibble, history)
#totalforecast
asdf |> group_by(month, .model) |> summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = .model)) +
geom_line() +
theme_minimal() +
labs(title = "Total") +
scale_color_manual(values = c("red","blue","black"))
asdf |>
group_by(month, .model, category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = .model)) +
geom_line() +
theme_minimal() +
labs(title = "Category") +
facet_wrap(vars(category), scales = "free", ncol = 1) +
scale_color_manual(values = c("red","blue","black"))
for (target_category in unique(data2$category)) {
p1 <- asdf |> filter(category == target_category) |>
group_by(month, .model, sub_category) |>
summarise(quantity = sum(quantity), .groups = "drop") |>
ggplot(aes(x = month, y = quantity, color = .model)) +
geom_line() +
theme_minimal() +
labs(title = target_category) +
facet_wrap(vars(sub_category), scales = "free") +
scale_color_manual(values = c("red","blue","black"))
print(p1)
}
6.5 Forecasting as integer values
When forecasting at the category level, we were okay rounding the forecast values to whole numbers, and dropping the few negative numbers to zero because we had high volume. The median forecast value for the new product_id level forecast is only .6. If we were just to round we would wipe out a lot of demand.
Depending on what the end use for the forecast is, we might need to fight a way to get these to whole numbers. I haven’t seen a case where MRP can accept anything else, so we’ll move to roll these up
item_fcst <- ttl_ag_wdrops_fit_rec_fc |>
as_tibble() |>
filter(!is_aggregated(product_id), .model == "ols") |>
mutate(category = as.character(category), sub_category = as.character(sub_category), product_id = as.character(product_id)) |>
select(month, category, sub_category, product_id, .mean)
x <- round(median(item_fcst$.mean), 3)
item_fcst |> ggplot(aes(x = .mean)) +
geom_histogram(binwidth = 0.1, col = "white", fill = "grey4") +
geom_vline(xintercept = x, col = "red", linewidth = 1) +
theme_minimal() +
labs(
title = paste("Median fcst value is ", x, " per month"),
x = "forecast"
) +
scale_x_continuous(breaks = seq(0,100,by = 1))
One method is to sum the sequences of forecasts, and make a whole number record once the cumulative sum reaches a new threshold. For example, a forecast of 0.25, 0.25, 0.25, 0.25, 0.25 would create a cumulative sum of 0.25, 0.5, 0.75, 1.0, 1.25, and the forecast values would be 0, 0, 0, 1, 0.
The code below accomplishes this by creating a cumulative sum for each forecast, a rolling calculation taking the floor value (e.g. floor of 1.9999 would be 1, and 2.000001 would be 2. Then a third column checking the difference at lag 1 and recording when it reaches a new integer threshold.
item_fct_int <- item_fcst |>
arrange(product_id, month) |>
group_by(product_id) |>
mutate(
csum = cumsum(.mean),
rolling_int = floor(csum),
int_fcst = rolling_int - lag(rolling_int, n = 1),
int_fcst = ifelse(is.na(int_fcst), 0, int_fcst)
)
item_fct_int |>
mutate(.mean = round(.mean, 2),
csum = round(csum, 2)) |>
select(month, product_id, .mean, csum, rolling_int, int_fcst) |>
head(24) |>
kable(caption = "Example with two products \n
.mean = forecast value \n
csum = cummulative sum of forecast vlues \n
rolling_in = floor value of csum \n
int_fcst = our forcast value, that takes rolling_int difference at lag 1",
booktabs = T)| month | product_id | .mean | csum | rolling_int | int_fcst |
|---|---|---|---|---|---|
| 2018 Jan | FUR-BO-10000330 | 0.16 | 0.16 | 0 | 0 |
| 2018 Feb | FUR-BO-10000330 | 0.11 | 0.27 | 0 | 0 |
| 2018 Mar | FUR-BO-10000330 | 0.21 | 0.48 | 0 | 0 |
| 2018 Apr | FUR-BO-10000330 | 0.29 | 0.77 | 0 | 0 |
| 2018 May | FUR-BO-10000330 | 0.30 | 1.07 | 1 | 1 |
| 2018 Jun | FUR-BO-10000330 | 0.37 | 1.44 | 1 | 0 |
| 2018 Jul | FUR-BO-10000330 | 0.23 | 1.66 | 1 | 0 |
| 2018 Aug | FUR-BO-10000330 | 0.19 | 1.86 | 1 | 0 |
| 2018 Sep | FUR-BO-10000330 | 0.39 | 2.25 | 2 | 1 |
| 2018 Oct | FUR-BO-10000330 | 0.33 | 2.58 | 2 | 0 |
| 2018 Nov | FUR-BO-10000330 | 0.48 | 3.06 | 3 | 1 |
| 2018 Dec | FUR-BO-10000330 | 0.62 | 3.67 | 3 | 0 |
| 2018 Jan | FUR-BO-10000362 | 0.34 | 0.34 | 0 | 0 |
| 2018 Feb | FUR-BO-10000362 | 0.29 | 0.63 | 0 | 0 |
| 2018 Mar | FUR-BO-10000362 | 0.39 | 1.02 | 1 | 1 |
| 2018 Apr | FUR-BO-10000362 | 0.47 | 1.49 | 1 | 0 |
| 2018 May | FUR-BO-10000362 | 0.48 | 1.97 | 1 | 0 |
| 2018 Jun | FUR-BO-10000362 | 0.55 | 2.51 | 2 | 1 |
| 2018 Jul | FUR-BO-10000362 | 0.41 | 2.92 | 2 | 0 |
| 2018 Aug | FUR-BO-10000362 | 0.37 | 3.29 | 3 | 1 |
| 2018 Sep | FUR-BO-10000362 | 0.57 | 3.86 | 3 | 0 |
| 2018 Oct | FUR-BO-10000362 | 0.51 | 4.37 | 4 | 1 |
| 2018 Nov | FUR-BO-10000362 | 0.66 | 5.03 | 5 | 1 |
| 2018 Dec | FUR-BO-10000362 | 0.80 | 5.82 | 5 | 0 |
Wow! To my surprise, the rolling int method actually worked pretty good aligning with our higher level forecasts. Only problem is that of course we have ZERO forecast for the first month, because there is no difference avail for the first month because there is now prior period!
This may not be a problem for supply chain planning if they have a frozen forecast zone (since we wouldn’t be able to change the forecast for one month ahead.
item_fct_int |>
group_by(month, category) |>
summarise(
dbl_fcst = sum(.mean, na.rm = T),
int_fcst = sum(int_fcst, na.rm = T),
.groups = "drop"
) |>
pivot_longer(cols = c("dbl_fcst", "int_fcst"), names_to = "fcst", values_to = "quantity") |>
ggplot(aes(x = month, y = quantity)) +
geom_line(aes(color = fcst)) +
theme_minimal() +
facet_wrap(vars(category), scales = "free") +
theme(axis.text.x = element_text(vjust = .5, angle = 45))
Its good practice to review top products. Ideally there is a set process but for this we’ll just take a look at the top 10. Below is history, our reconciled hierarchical forecast as a decimal “dbl_fcst” and the rolling integer version we just created (int_fcst).
Our integer forecast is aligning pretty nicely with the decimal forecast, but the forecast might be going off the rails a little. This would be a good follow up opportunity with the sales team to figure out what the root cause of increased sales are.
Recall back to our histogram in the data exploration, these high monthly volumes in 2017 are way outside the normal sales and clearly unusual.
top9_list <- ttl_ag_wdrops_fit_rec_fc |> as_tibble() |>
filter(!is_aggregated(product_id), .model == "ols") |>
select(product_id, .mean) |>
mutate(product_id = as.character(product_id)) |>
group_by(product_id) |>
summarize(quantity = sum(.mean), .groups = "drop") |>
arrange(desc(quantity)) |>
top_n(9, quantity)
a <- ttl_ag_wdrops_fit_rec_fc |> filter(!is_aggregated(product_id), .model == "ols") |>
as_tibble() |>
mutate(product_id = as.character(product_id),
quantity = .mean) |>
filter(product_id %in% top9_list$product_id) |>
select(product_id, month, .model, quantity) |>
rename(type = .model) |>
mutate(type = "dbl_fcst")
b <- item_fct_int |>
filter(product_id %in% top9_list$product_id) |>
mutate(type = "int_fcst") |>
select(product_id, month, type, int_fcst) |>
rename(quantity = int_fcst)
c <- data2_ts_monthly |> as_tibble() |>
mutate(type = "history") |>
filter(product_id %in% top9_list$product_id) |>
select(product_id, month, type, quantity)
asdf2 <- rbind(a,b, c)
for (target_prod in unique(asdf2$product_id)) {
p1 <- asdf2 |>
filter(product_id == target_prod) |>
ggplot(aes(x = month, y = quantity, color = type)) +
geom_line() +
theme_minimal() +
labs(title = paste0(target_prod, " Forecast")) +
scale_color_manual(values = c("red", "black", "blue"))
print(p1)
}
